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Abstract 

Probiotics are live microorganisms that potentially confer beneficial outcomes to host by modulating 
gut microbiota in the intestine. The aim of this study was to comprehensively investigate effects of probio- 
tics on human intestinal microbiota using 454 pyrosequencing of bacterial 1 6S ribosomal RNA genes with 
an improved quantitative accuracy for evaluation of the bacterial composition. We obtained 1 58 faecal 
samples from 1 8 healthy adult Japanese who were subjected to intervention with 6 commercially available 
probiotics containing either Bifidobacterium or Lactobacillus strains. We then analysed and compared bac- 
terial composition of the faecal samples collected before, during, and after probiotic intervention by 
Operational taxonomic units (OTUs) and UniFrac distances. The results showed no significant changes 
in the overall structure of gut microbiota in the samples with and without probiotic administration regard- 
less of groups and types of the probiotics used. We noticed that 32 OTUs (2.7% of all analysed OTUs) 
assigned to the indigenous species showed a significant increase or decrease of > 1 0-fold or a quantity dif- 
ference in >1 50 reads on probiotic administration. Such OTUs were found to be individual specific and 
tend to be unevenly distributed in the subjects. These data, thus, suggest robustness of the gut microbiota 
composition in healthy adults on probiotic administration. 
Key words: probiotics; gut microbiota; 1 6S ribosomal RNA gene; pyrosequencing 



1 . Introduction 

Probiotics are defined as live bacterial strains con- 
ferring various benefits to the consumer by modulat- 
ing the intestinal ecosystem, thereby potentially 
promoting host health and improving host disease 
risk. 1-11 Various probiotic strains have been industri- 
ally developed and marketed as a variety of products 
and applications such as fermented foods and 



supplements, including yogurt 12-15 Most probiotics 
taxonomically belong to two genera, Bifidobacterium 
and Lactobacillus, that originate from various environ- 
ments, including the human intestine, and both 
species are generally regarded as safe. 16-18 

The interaction between administrated probiotics 
and indigenous microbiota is one of the most attract- 
ive and important research areas, particularly because 
gut microbiota have been shown to be profoundly 
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associated with various host physiology states, includ- 
ing disease, diet, and age through the shift of bacterial 
composition, as well as metabolic and nutritional pro- 
cesses. 19-23 The ability of probiotics to survive 
through the intestine and to modulate gut microbiota 
is a critical factor in determining their potential for 
health-related outcomes. 

There have been a large number of probiotic inter- 
vention studies to assess the impact of probiotics on 
gut microbiota in healthy adults, 24-34 infants, and 
children, 35,36 and in clinical trials on patients with a 
variety of diseases. 37,38 Most probiotic intervention 
studies were carried out by comparison between pro- 
biotic-treated groups and placebo controls and exam- 
ined only one or two samples from periods before and 
during intervention or post-intervention for each 
subject. These experimental designs make evaluation 
of results obscure from a statistical viewpoint due to 
the high inter-individual variability of gut micro- 
biota. 4 In addition, most of the analyses focussed on 
the composition of specific bacterial species or 
groups by conventional methods such as culturing, 
quantitative polymerase chain reaction (qPCR), fluor- 
escence in situ hybridization, denaturing gradient gel 
electrophoresis, or terminal-restriction fragment 
length polymorphism based on the bacterial 1 6S 
ribosomal RNA gene (1 6S). These conventional meth- 
odologies may also overlook subtle changes in bacter- 
ial community structure and change of species other 
than targeted species. Thus, the effect of probiotic ad- 
ministration on the overall structure of gut microbiota 
is largely unknown. 

Recently, a high-throughput sequencing-based ana- 
lysis has been conducted for gut microbiota fed with a 
probiotic yogurt that provided new insights into pro- 
biotics research by utilizing a large-scale dataset. 39 
However, much more data are required to understand 
the impact of probiotics on gut microbiota. Recent 
advances in sequencing technology have enabled us 
to elucidate complex bacterial communities, includ- 
ing human gut microbiota. 40,41 Particularly, 454 pyr- 
osequencing of bacterial 1 6S gene tags coupled with 
bioinformatics provides a high-throughput and cost- 
effective approach for the comprehensive analysis of 
bacterial communities at the species level. 42-48 

In this study, we developed an analysis pipeline for 
bacterial communities based on barcoded 454 pyro- 
sequencing of 1 6S gene tags using modified PCR 
primers that improved the quantitative accuracy of in- 
ferred species composition in human gut microbiota. 
Using this pipeline, we analysed faecal samples longi- 
tudinally collected from individuals with and without 
probiotic administration to evaluate the effect of 
probiotics on gut microbiota with respect to species 
richness and diversity. The results revealed the 



robustness and stability of gut microbiota of healthy 
adults in response to probiotic administration. 

2. Materials and methods 

2.1 . Subjects, faecal sample collection, and probiotic 
intervention 

Eighteen healthy volunteers (age: 22 + 3.1 6 yrs, 6 
male, 1 2 female) were recruited through Azabu 
University, Kanagawa, Japan (Supplementary Table S1). 
All subjects were informed of the purpose of this study. 
This study was approved by the ethical committee of 
Azabu University, and written consent was obtained 
from all subjects. No subjects were treated with antibio- 
tics during faecal sample collection. The subjects were 
divided into six groups (three subjects per group), and 
each group consumed six different commercially avail- 
able probiotics supplied from Yakult Honsha Co., Ltd, 
Kagome Co., Ltd, Morinaga Milk Industry Co., Ltd, 
Takanashi Milk Products Co., Ltd, Meiji Co., Ltd, and 
Danone Japan Co., Ltd, respectively (Supplementary 
TableS! ).The numberof each bacterial strain contained 
in the probiotic products was estimated as the genome 
equivalent by qPCR of 1 6S ribosomal RNA genes using 
2 7Fmod-33 8R, followed by pyrosequencing of the 
1 6S amplicons (see below). The genome equivalent 
per gram or millilitre and the total genome equivalent 
of each bacterial strain in one probiotic product are sum- 
marized in Supplementary Table S1. Three subjects in 
each group consumed the same probiotics daily for 
8 weeks according to the schedule of sampling and pro- 
biotic intervention (Supplementary Fig. S1). Faecal 
samples from 4 weeks before (S00) and 8 weeks 
during probiotic intervention (S01 -S04), and 8 weeks 
after cessation of probiotic intervention (S05-S08), 
were collected every 2 weeks from each subject. In 
total, we collected 1 58 faecal samples from the 1 8 sub- 
jects because we could not collect 1 sample each from 4 
of the subjects. 

2.2. Recovery of bacteria from faecal samples 
Freshly collected faeces (1 .0 g) were suspended in 

20% glycerol (Wako Pure Chemical Industries, Ltd) 
and phosphate buffered saline solution (Life 
Technologies Japan, Ltd, Tokyo, Japan), frozen in 
liquid nitrogen, and stored at -80°C until ready for 
use. Bacterial pellets were prepared from frozen 
faecal samples as described previously. 49 

2.3. DNA isolation from bacteria 

Faecal DNA was isolated and purified according to 
the literature, with minor modifications. 49 The bac- 
terial pellet was suspended and incubated with 
1 5 mg/ml lysozyme (Sigma-Aldrich Co., LCC) at 
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37°C for 1 h in TE1 0. Purified achromopeptidase 
(Wako Pure Chemical Industries, Ltd) was added at a 
final concentration of 2000 units/ml and then incu- 
bated at 37°C for 30 min. The suspension was 
treated with 1% (wt/vol) sodium dodecyl sulphate 
and 1 mg/ml proteinase K (Merck Japan) and 
incubated at 55°C for 1 h. The lysate was treated 
with phenol/chloroform/isoamyl alcohol (Life 
Technologies Japan, Ltd). DNA was precipitated by 
adding ethanol and pelleted by centrifugation at 
3,300 g at 4°C for 1 5 min. The DNA pellet was 
rinsed with 75% ethanol, dried, and dissolved in 
10 mM Tris-HCI/1 mM EDTA (TE). DNA samples 
were purified by treating with 1 mg/ml RNase A 
(Wako Pure Chemical Industries, Ltd) at 37°C for 
30 min and precipitated by adding equal volumes of 
20% polyethylene glycol solution (PEG6000-2.5M 
NaCl). DNA was pelleted by centrifugation at 8,060 
g at 4°C, rinsed with 7 5% ethanol, and dissolved in TE. 

2.4. 454 barcoded pyrosequencing of 1 6S rRNA gene 
The V1 -V2 region of the 1 6S rRNA gene was amp- 
lified using forward primer (5'-CCATCTCATCCCTGCG 
TGTCTCCGACTCAG NNNNNNNN N Nagrgtttgaty mtggc 
tcag-3') containing the 454 primer A, a unique 1 0-bp 
barcode sequence for each sample (indicated in N), 
and 2 7Fmod (5'-agrgtttgatymtggctcag) in which the 
third base A in the original primer 27F was changed 
to R, and reverse primer (5'-CCTATCCCCTGTGTGC 
CTTGGCAGTCTCAGtgctgcctcccgtaggagt-3') contain- 
ing the 454 primer B and reverse primer 338R (5'- 
tgctgcctcccgtaggagt). PCR was performed in 1 x Ex 
Taq PCR buffer (50 jjiI), deoxynucleoside triphosphate 
(2.5 mM), Ex Taq polymerase (Takara Bio, Inc., Shiga), 
each primer (1 0 |xM), and 40 ng of extracted DNA 
under conditions of 2 min at 96°C, 20 cycles of 
96°C for 30 s, 55°C for 45 s, and 72°C for 1 min, 
and a final extension of 72°C for 1 0 min on a 9700 
PCR system (Life Technologies Japan, Ltd, Tokyo, 
Japan). PCR products of approximately 370 bp were 
confirmed by agarose gel electrophoresis, purified by 
AMPure XP magnetic purification beads (Beckman 
Coulter, Inc., Brea, CA, USA), and quantified using 
the Quant-iT PicoGreen dsDNA Assay Kit (Life 
Technologies Japan, Ltd, Tokyo, Japan). Mixed 
samples were prepared by pooling approximately 
equal amounts of PCR amplicons from each sample 
and subjected to 454 GS FLX Titanium or 454 GS 
JUNIOR (Roche Applied Science) sequencing accord- 
ing to the manufacturer's instructions. 

2.5. Analysis pipeline for 454 barcoded 
pyrosequencing of 1 6S PCR amplicons 

We developed an analysis pipeline for 454 bar- 
coded pyrosequencing of PCR amplicons of the V1 -2 



region amplified by 27Fmod-338R primers. First, 

1 6S reads were assigned to each sample based on 
the barcode sequence information. Second, 1 6S 
reads that did not have PCR primer sequences at 
both sequence termini and those with an average 
quality value < 25 were filtered out. Third, 1 6S 
reads containing possible chimaeric sequences that 
had BLAST match lengths of<90% with reference 
sequences in the database were removed. Reads 
removed in these processes accounted for about 
35% of all reads, most of which represented reads 
lacking PCR primer sequences (Supplementary Table 
S2). Finally, filter-passed reads were obtained for 
further analysis by trimming off both primer 
sequences. 

All 3000 filter-passed reads of the 1 6S V1-2 
sequences obtained from each subject were deposited 
in DDBJ/GenBank/EMBL with accession numbers 
DRA000869-DRA000886. 

2.6. Assessment of the quantitative accuracy of 1 6S 
data using artificial bacterial communities 

Two artificial bacterial communities (designated 
'mockOV and 'mocl<02') were constructed by mixing 
genomic DNA from 1 0 and 1 1 different human gut- 
associated bacterial strains with an appropriate ratio, 
respectively (Supplementary Table S3). Genome 
sequences of these microbes were completely 
sequenced and are publicly available. From these 
communities, we amplified the V1 -2 region by PCR 
using 27F-338R and 27Fmod-338R primers, the 
V5-6 region by 787F-1061R primers, and the V1-9 
region by 27F-1 492R primers. V1 -2 and V5-6 ampli- 
cons were subjected to 454 pyrosequencing, and V1- 
9 amplicons were cloned in Escherichia coli, and 3000 
clones were sequenced by the Sanger method, and 
the products were analysed with the ABI3730xl (Life 
Technologies Japan, Ltd, Tokyo). We also performed 
duplicate qPCR experiments targeting a specific 
genomic region of the bacterial strains in the two 
mock communities. All filter-passed 1 6S de novo 
sequences and qPCR data were then analysed by prin- 
cipal component analysis (PCA) to compare and 
assess the quantitative accuracy (Fig. 1 ). 

The error rate of the filter-passed sequences using 

2 7Fmod-338R primers obtained from the two mock 
communities was estimated by aligning the 1 6S V1- 
2 de novo sequences with the reference 1 6S sequences 
in the two mock communities (Supplementary 
Table S4). 

2.7 . Data analysis 

2.7.1. Database Two databases were con- 
structed for the analysis of 1 6S sequences. One is 
the 1 6S rRNA gene sequence database constructed 
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Figure 1. Assessment of the quantitative accuracy of the analysis of the bacterial composition of two mock communities by various 
methods. PCA analysis of the data was obtained from various methods using rnockOI (a) and mock02 (b). Closed circle: expected, 
open circle: duplicate qPCR, closed square: pyrosequencing of 1 6S V1-2 region using 27Fmod, open square: pyrosequencing of 1 6S 
V1-2 region using 27F, closed triangle: pyrosequencing of 1 6S V5-6 region, open diamond: Sanger sequencing of nearly full-length 
1 6S clone. 



by collecting 1 6S sequences of >1 200 bp of bacterial 
isolates from the Ribosomal Database Project 
v. 1 0.2 7. Another database is the reference genome 
database constructed by collecting genome 
sequences from the NCBI FTP site (ftp://ftp.ncbi.nih. 
gov/genbank/, Dec 201 1) that includes 1482 com- 
plete and 605 draft bacterial genomes. 

2.7.2. Operational taxonomic unit (OTU) and 
UniFrac distance analysis We used 3000 

filter-passed reads of 1 6S sequences for operational 
taxonomic unit (OTU) and UniFrac distance analysis 
for each sample. For OTU analysis, clustering of 1 6S 
reads was done using a 96% pairwise-identity cutoff 
with the UCLUST program (www.drive5.com). 
Representative sequences for each OTU were assigned 
to bacterial species by BLAST search with a 96% pair- 
wise-identity cutoff against the two databases men- 
tioned above. UniFrac distance analysis was used to 
determine the dissimilarity (distance) between two 
communities based on the fraction of branch length 
shared between two communities within a phylogen- 
etic tree constructed from 1 6S sequence datasets. 44 

2.7.3. Other Estimation of OTU numbers by ex- 
trapolation (Chao! and ACE) was calculated with 
the vegan package (v2.0-5) for R (v2.1 5.2). 



3. Results and discussion 

3.1 . Quantitative accuracy of 1 6S data produced 
by 454 pyrosequencing 
Pyrosequencing of PCR amplicons of bacterial 1 6S 
short variable regions is the most popular and a 



high-throughput approach to infer and characterize 
the species composition in bacterial communi- 
ties. 42,45,46,48 The 454 pyrosequencing platform, 
which can produce over 400 bases per read, is also su- 
perior to shorter read-length sequencers with respect 
to sequence accuracy for single-end sequencing. 50,51 
However, this PCR-based method has a problem par- 
ticularly in quantification of the composition of the 
genus Bifidobacterium, a dominant species in human 
gut microbiota because the 1 6S sequence of 
Bifidobacterium has a few base mismatches with the 
commonly used PCR primer 2 7F (or 8F), underesti- 
mating this genus in the community. 52-55 To 
improve this, we modified primer 2 7F to 2 7Fmod 
by changing the third base G to R (G or A) in 2 7F- 
YM 53 that perfectly matched with the annealing site 
of the Bifidobacterium 1 6S gene (see Materials and 
methods). 

To assess the 1 6Sdata using 2 7Fmod, we compared 
various 1 6S sequence and qPCR data obtained from 
two mock communities (Supplementary Table S3) 
that are useful to evaluate the quantitative accuracy 
of 1 6S-based data and the sequencing error 
rate. 56,57 Quantitative accuracy of the overall bacterial 
composition was evaluated by comparing the similar- 
ity of each data to the expected ('Expected') using PCA 
(Fig. 1). From the PCA data, Euclidean distance was 
calculated for evaluation of the similarity of each 
data with the 'Expected'. The results revealed that 
the order of their similarities with the 'Expected' was 
the qPCR data > the V1 -2 data using 2 7Fmod > the 
V5-6 data > the V1-2 data using 27F>>the data 
of Sanger sequencing-based full-length V1 -9, indicat- 
ing that the use of 27Fmod greatly improved the 
quantitative accuracy for evaluation of the overall 
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bacterial composition (Supplementary Table S5). This 
improvement was largely dependent on the improved 
estimation of the Bifidobacterium content by the use 
of 2 7Fmod. The average relative Bifidobacterium 
content in the two mock communities estimated 
from the data of V1-2 using 2 7F was only 1.5% of 
the 'Expected' (100%), whereas the use of 27Fmod 
increased the relative content to 61% that was also 
better than that estimated from the data of V5-6 
and Sanger full-length analyses (Supplementary Fig. 
S2). Because qPCR can be used only when genomes 
of all bacteria in a given community are known, or 
only for a limited number of specific known species, 
we concluded that 454 pyrosequencing of the V1-2 
region using 27Fmod-338R provided more quantita- 
tively accurate data for bacterial composition in 
human gut microbiota than that using the conven- 
tional 2 7F primer. 

We estimated the average error rate of filter-passed 
V1 -2 data using 27Fmod-338R by aligning the V1 -2 
and reference 1 6S sequences of bacterial strains 
used in the two mock communities. The error rate 
was estimated to be 0.58 and 0.66% for mockOI 
and mocl<02 by local alignment, respectively 
(Supplementary Table S4). These error rates are 
similar to the previously published data, 43,45,50 but 
lower than in another study. 58 The latter may be 
due to differences in the examined alignment length 
and between local and global alignments. Errors in 
454 pyrosequencing data can be the primary cause 
for overestimation of the OTU number that is an 
issue which needs to be improved for accurate estima- 
tion of species richness in bacterial communities. 59,60 
We compared OTU numbers generated from cluster- 
ing of various qualities of 1 6S reads with a 96% and 
a 97% pair-wise identity cutoff. For this comparison, 
we made and used three datasets: only primer 
check-passed reads having the highest error rates, 
filter-passed reads, and selected filter-passed reads 
having the lowest error rates. The results indicated 
that a 96% cutoff clustering of error-rich reads and 
a 97% cutoff clustering of filter-passed reads gave 
the worse results than a 96% cutoff clustering of 
filter-passed and selected filter-passed reads 
(Supplementary Fig. S3). A 97% cutoff was defined 
for clustering of highly accurate Sanger full-length 
1 6S sequences. 61 Therefore, in clustering of pyrose- 
quencing data having higher error rate than Sanger 
data, the use of a cutoff identity lower than 97% 
and a lower number of reads are reasonable to 
reduce overestimation of the OTU number. A 96% 
cutoff clustering of filter-passed reads gave similar 
OTU numbers up to 30-50 reads to those of filter- 
passed reads having the lowest error rates. These 
read numbers are approximately three to five times 
the number of input strains. After several trials 



testing the mock communities, we decided to use 
3000-5000 reads per sample for clustering with a 
96% cutoff for the analysis of human gut microbiota. 
Indeed, OTU numbers using a 96% cutoff clustering of 
3000 reads decreased about 15% when compared 
with those using a 97% cutoff clustering. 

3.2. Species richness and diversity in human faecal 
microbiota with probiotic intervention 
We randomly selected 3000 reads of 1 6S V1-2 
sequences from all filter-passed reads for each 
sample (Supplementary Table S2) and used 
474 000 reads in total from 158 faecal DNA 
samples of 1 8 subjects for the analysis of species rich- 
ness and composition in human gut microbiota. 
Clustering of all reads with a 96% pairwise-identity 
cutoff gave a total of 2758 OTUs. After removing 
the minority OTUs having <0.1% abundance in any 
samples, 1 175 OTUs having >0.1% abundance in at 
least 1 sample, accounting for 99.1% of all 1 6S 
reads, were used for further analysis. 

3.2.1 . Detection of administrated probiotic strains in 
faecal samples We investigated whether 
administrated strains contained in the probiotic pro- 
ducts can be detected in faecal DNA. We sequenced 
the 1 6S V1 -2 region of all bacterial strains contained 
in probiotic products used in this study. The BLAST 
search to the databases indicated that except for the 
Bifidobacterium long urn strain used in Group III, the 
1 6S sequences of all strains in the probiotic products 
significantly differed from those of the indigenous 
species phylogenetically closest to the probiotic 
strains. The 1 6S sequence of the B. longum strain 
used in Group III was almost identical to that of an in- 
digenous Bifidobacterium species, so that we used a 
distinguishable additive Lactococcus lactis strain in 
this product for the detection of administrated bac- 
teria in Group III samples. The 1 6S sequences of 
these probiotic strains were included in the databases 
constructed in this study, and the 1 6S reads assigned 
to administrated strains had the average similarity 
between 99.4 and 99.9% identities with the reference 
sequences (data not shown). The 1 6S reads assigned 
to probiotic or additive strains were detected in 
samples (S01-S04) during probiotic intervention 
[designated Tro(+)'] at various frequencies, but 
almost none were detected in samples (S00 and 
S05-S08) without probiotic administration [desig- 
nated 'Pro(-)'] (Supplementary Table S6). The admi- 
nistrated probiotic strains were shown to be more 
frequently detected in samples during the interven- 
tion than in the pre- and post-intervention periods 
using different detection methods such as culturing, 
targeted PCR, and hybridization. 24,26 " 28,30,32 ' 33,62 In 
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the present study, two probiotic Lactobacillus and one 
additive Lactococcus strains were detected in post- 
intervention samples in three subjects with a 
minimum count, respectively. The similarity of three 
1 6S sequences was 99.4, 99.7, and 100% identity 
with those of administrated Lactobacillus and 
Lactococcus strains, indicating that these are admini- 
strated strains. The survival of some probiotic strains 
in the post-intervention period was also reported 
previously. 28,30 Our data suggested that some pro- 
biotic strains seem to be able to persistently colonize 
the intestine and their survivability may be related to 
metabolic activity in the intestine. 63,64 Probiotic 
Bifidobacterium strains were not detected in any 
Pro(-) samples. However, we found two distinct 1 6S 
sequences both assigned to Bifidobacterium animalis 
in two subjects APr37 and APr39. One showed a 
high similarity of >98% identity with the 1 6S se- 
quence of the administrated B. animalis and was 
detected with high frequency only in the Pro(+) 
samples, whereas another showed a low similarity of 
96.5-97.4% identity (a mean of 97.2%) with low fre- 
quency in both the Pro(-) and Pro(+) samples. These 
data suggest the presence of unknown indigenous 
species phylogenetically close to, but distinct from, 
probiotic B. animalis in human gut microbiota. The 
total number of bacteria contained in each probiotic 
product was varied between 1 0 9 and 1 0 10 , showing 
no large difference in quantity among them 
(Supplementary Table S1). No clear correlation was 



also observed between the number of bacteria in 
the products and the frequency in detection of the 
administrated strains in the Pro(+) samples. From 
these observations, the frequency of administrated 
bacteria detected in faeces may not be largely affected 
by their amounts in the products. Therefore, detec- 
tion of Lactobacillus brevis and Lactobacillus delbrueck 
at relatively low level in faeces cannot be simply 
explained by the difference in a dose, but could be 
considered the association with several factors such 
as their survivability in the intestine, diet, or physio- 
logical conditions of subjects. 

3.2.2. Change of species richness in samples with 
and without probiotics We analysed 
species richness (OTU number) in the Pro(+) and 
Pro(-) samples. Supplementary Figure S4 shows the 
change in OTU numbers for every sample in each 
subject, indicating that OTU numbers vary dramatic- 
ally for every sample. Most of the variation can be 
attributed to single OTUs representing the minority 
species. We averaged the OTU numbers of the 
Pro(-) and Pro(+) samples and compared them for 
subject, group, type of probotics (Lactobacillus and 
Bifidobacterium), and all combined samples, respect- 
ively (Fig. 2). The average OTU numbers in 6 out of 
1 8 subjects were decreased in the range of the ratio 
of 0.83-0.95 in the Pro(+) samples when compared 
with the Pro(-) samples, whereas those in other 1 2 
subjects were increased in the range of the ratio of 




Figure 2. Change in OTU number in faecal microbiota with and without probiotic administration, (a) Individual, (b) group, (c) type of 
probiotics. Black bar indicates Pro(-) samples. Grey bar indicates Pro(+) samples. The error bars represent standard deviation. 
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1.01-1.43. For group, only Group IV showed a de- 
crease in the average OTU number in the Pro(+) 
samples with the ratio of 0.94. For type of probiotics 
and all samples, the average OTU numbers in the 
Pro(+) samples were slightly more abundant (ap- 
proximately 1.07-fold) than those in the Pro(-) 
samples, but no statistical significance was observed 
in any dataset. The increase in OTU number in the 
Pro(+) samples was largely due to the minority 
species (Supplementary Fig. S4), whereas the abun- 
dance of the majority species (OTUs containing >1 0 
reads) was almost constant over time. We performed 
the same analysis using different sets of 3000 reads 
for each subject. The analysis reproducibly showed 
the similar pattern and the degree of the change in 
OTU numbers to which the minority species is 
largely attributed (data not shown). These data indi- 
cate that administration of probiotics tends to in- 
crease species richness in faecal microbiota that may 



be beneficial for the consumer because the species 
richness in faecal microbiota of subjects afflicted 
with disease such as inflammatory bowel disease is 
significantly reduced when compared with that of 
healthy subjects. 65 

3.2.3. Change of species composition in samples 
with and without probiotics We obtained 
the average weighted and unweighted UniFrac dis- 
tances within Pro(-), within Pro(+), and between 
Pro(-) and Pro(+) samples for every group, probiotic 
types, and all subjects, respectively (Fig. 3). High 
UniFrac distance implies high variability of microbiota 
structure within and between samples. If the differ- 
ence between any pair of the three distances is statis- 
tically significant, it can be considered that probiotic 
administration significantly affected the overall micro- 
biota composition. We found the largest difference 
between weighted UniFac distances of the Pro(+) 
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Figure 3. Average UniFrac distance within Pro (-) and Pro(+) and between Pro(-) and Pro(+) for each group, type of probitotics, and all 
subjects. Average UniFrac distance between any pair of the three distances for type of probiotics and all subjects (a and b), and each 
group (c and d). Open circle, closed circle, and closed square indicate average UniFrac distance within Pro (-), within Pro (+), and 
between Pro(-) and Pro(+) samples, respectively. 
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and Pro(-) samples in Group VI. However, statistical 
evaluation of this difference by the Student's f-test 
showed no significance (P-value > 0.05) for 781 out 
of 1 000 times (Supplementary Table S7). These 
data imply high stability of gut microbiota to probiotic 
administration for all subjects examined. We also ana- 
lysed UniFrac distances of intra-subject gut microbiota 
(Fig. 4). Although 5 subjects (APr02, 1 2, 1 6, 3 7, and 
39) showed a significant difference in the UniFrac 
distances between Pro(-) and Pro(+) samples, the 
resultsshowed that both weighted and unweighted dis- 
tances between Pro(-) and Pro(+) of all intra-subjects 
were significantly lower than the average distance 
of the 18 unrelated subjects. The Welch's Mest 
for these differences showed statistical significance 
(Supplementary Table S8). We also performed the 
UniFrac distance analysis using different 1 6S datasets 
of 5000 reads for group, type of probiotics, all subjects, 
and intra-subject. The results similarly showed no stat- 
istical significance in differences between any pair of 
the 3 UniFrac distances and the significantly lower 
UniFrac distance of each intra-subject than that of the 
1 8 unrelated subjects (data not shown). Thus, these 
data suggested that the perturbation of microbiota eli- 
cited by probiotics in an intra-subject did not overcome 
the inter-subject variations of gut microbiota, support- 
ing high intra-specificity and stability of gut micro- 
biota. 66,67 This robustness of gut microbiota of adults 
is in contrast with the profound effect of antibiotic 



administration on adult gut microbiota 68 and the 
observed response of gut microbiota of infants fed 
with probiotics, in which the infant gut microbiota 
composition was considerably affected by probiotics. 36 
A short-term dietary intervention study showed that in 
controlled feeding of the same diet to subjects over 1 0 
days, a marked change was observed within 1 day after 
the intervention initiation. 69 In the present study, no 
significant difference was observed between samples 
before (S00) and first samples (S01 ) afterthe interven- 
tion initiation (data not shown). It would be valuable to 
analyse faecal samples collected within a few days after 
administration of probiotics for evaluation of the short- 
term effect of probiotics. 



4. Identification of bacterial species showing 
significant increase or decrease by probiotic 
administration 

Although our results suggested that administration 
of probiotics had almost no effect on the overall 
structure of gut microbiota, it is possible to identify 
bacterial species largely responding to the admini- 
strated probiotics at the OTU/species level. We sur- 
veyed OTUs showing an increase or a decrease 
between the Pro(+) and Pro(-) samples by comparing 
the number of 1 6S reads for each OTU. We first enum- 
erated the OTUs showing > 2-fold change between the 



u 
u 


0.6 n 


c 




a 


0.5- 








0.4- 






H 


0.3- 




0.2 


■a 






s 

M 


0.1 - 


W( 


0 



(a) 



O 



Weighted 



I 



O 



I 

o 



■ 

o 



o 



<*y *y <*r 

V V V V 



<6 



or) »■ 
*?• *5- 4& 



0» 

g 0.8 

1 0.7 
•- 

= 0.6 



0.5 



0.4 



(b) 

B 



o 
■ 



Unweighted 



! a . 

•of 



(VJ v v 



o 



■ 

o 



I 

o 



I 1 



e 



I • 

o 



o o 
■ ■ 



Figure 4. Average UniFrac distance within Pro(-) and Pro(+) and between Pro(-) and Pro(+) for each subject. Open circles, closed circles, 
and closed squares indicate average UniFrac distance within Pro(-), within Pro(+), and between Pro(-) and Pro(+) samples, 
respectively. Closed triangles indicate average UniFrac distance between samples (S00) of 1 8 unrelated individuals. 
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Pro(-) and Pro(+) samples for each subject, and the 
quantity difference was also obtained by subtracting 
the 1 6S read number of the Pro(+) samples from 
that of the Pro(-) samples. This is because OTUs 
showing a high quantity difference, but less fold 
change may also have substantial influence on gut 
microbiota composition. We found several OTUs sig- 
nificantly changed by probiotic administration, includ- 
ing OTUs assigned to both the indigenous and 
administrated strains (Fig. 5). We listed 88 OTUs 
(7.5% of all analysed 1175 OTUs) showing significant 
change of >3-fold, among which 30 OTUs changed 
by >1 0-fold (Supplementary Fig. S5). We excluded 6 
OTUs assigned to the administrated strains from the 
30 OTUs and obtained 24 OTUs assigned to the indi- 
genous species, including OTU00072 assigned to 
Streptococcus salivarius that showed significant 
change in 2 subjects (Supplementary Table S9). We 
also found seven OTUs showing significant difference 
in quantity between both samples (Supplementary 
Table S1 0). Of the combined 32 OTUs (2.7%), 18 
were increased and 1 4 were decreased by probiotic ad- 
ministration. Many of the OTUs showing a significant 
increase were assigned to minority species in the 
Pro(-) samples, but some increased up to nearly 7% 
in abundance (e.g. OTU00372 assigned to 



Eubacterium rectale). On the other hand, the OTUs 
showinga significant decrease were almost undetected 
in the Pro(+) samples. Phylum-level species 
assignment showed that species belonging to the 
phylum Firmicutes were most largely affected by both 
probiotics, and all species belonging to the phylum 
Bacteroidetes were affected only by Lactobacillus 
probiotics (Table 1). The 32 OTUs were assigned to 
2 7 indigenous species, among which 4 species 
(Clostridium clostridioforme, Eubacterium eligens, 
E. rectale, and Faecalibacterium prausnitzii) were 
assigned by 8 different OTUs and 1 species (S. salivarius) 
was assigned by the 2 same OTUs as described above. 
All these species except for S. salivarius were found to 
show significant change only in one subject, indicating 
that response of the indigenous species to probiotics is 
highly individual specific (Supplementary Fig. S6). 
Two different OTUs (OTU02677 and OTU02748) 
assigned to F. prausnitzii, of which the reduction is 
known to be correlated with inflammatory bowel 
disease, 70 were found to both decrease and increase 
in the same subject (APr40) by probiotic administra- 
tion, suggesting that these two phylogenetically close 
species may have the diversity of response to probiotic 
action. We also examined distribution of the 3 2 OTUs in 
the subjects. The results revealed that 4 subjects 
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Table 1. Phylum-level species assignment of OTUs showing significant fold change or quantity difference by administration of probiotics 



Type of 
probiotics 


Change 


a Number of 
varied OTUs 




Fold change (>1 0-fold) 




Number of 
varied 
w 1 US 


Difference (> 1 50 reads) 


Firmicutes 


Actinobacteria 


Bacteroidetes 


Unclassified 
bacterium 


Firmicutes 


Bacteroidetes 


Lactobacillus 


Increase 


9 


7 


0 


1 


1 


3 


2 


1 




Decrease 


7 


3 


1 


3 


0 


1 


1 


0 




Total 


1 6 


1 0 


1 


4 


1 


4 


3 


1 


Bifidobacterium 


Increase 


5 


5 


0 


0 


0 


1 


1 


0 




Decrease 


4 


4 


0 


0 


0 


2 


2 


0 




Total 


9 


9 


0 


0 


0 


3 


3 


0 


All 


Increase 


1 4 


1 2 


0 


1 


1 


4 


3 


1 




Decrease 


1 1 


7 


1 


3 


0 


3 


3 


0 




Total 


25 


19 


1 


4 


1 


7 


6 


1 



Administrated probiotic strains were excluded, and only OTUs with a P-value < 0.05 are shown. 
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increased and decreased OTUs, respectively. 



(APr1 1,2 0, 22, and 3 7) did not have such OTUs and 8 
subjects had only 1 OTU, whereas 4 subjects (APr01 , 
02, 31, and 39) had more than 4 OTUs showing signifi- 
cant change (Fig. 6), suggesting their uneven distribu- 
tion in the 1 8 subjects. These data imply existence of 
the sensitive and less sensitive responders to probiotic 
action and if so, it would be interesting to investigate 
the relation between gut microbiota type and its re- 
sponse to probiotics. 

In summary, we analysed changes of the gut micro- 
biota composition of healthy adults fed with probio- 
tics using the 454 pyrosequencing platform with the 
improved quantitative accuracy for evaluation of the 
overall bacterial composition. The present study 
using large datasets enabled us to more comprehen- 
sively and precisely evaluate the effect of probiotics 
on gut microbiota than the previous probiotic inter- 
vention researches in which the analysis exclusively 



focussed on only several limited bacterial species 
using conventional methods. Our data further 
support the high inter-subject variability and the 
high intra-subject stability that is the current 
common view for the feature of adult gut microbiota. 
A recent study of gut microbiota in twins demon- 
strated that probiotics had almost no effect on the 
community structure, but affected the gene expres- 
sion of microbiota. 39 To more deeply understand 
the potential function of probiotics, the analysis of 
bacterial and host cell's transcriptome and intestinal 
metabolome is required. 
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